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The numerical evaluation of slowly convergent integrals with infinite 
limits and regularly oscillating analytic integrands is discussed. Functions 
are presented that can be subtracted from the integrands to increase the 
rate of convergence. One of the integrals used as an example is an often- 
studied Fourier integral related to the power spectrum of a phase-modu- 
lated wave. It is pointed out that, when Fourier integrals whose values are 
known to be positive, e.g., probability densities, are evaluated by the 
trapezoidal rule, the values given by the trapezoidal rule are never less than 
the true values. 

I. INTRODUCTION 

This paper is in the nature of an addendum to an earlier paper in 
this journal dealing with the numerical evaluation of definite integrals. 1 

Integrals with infinite limits and oscillating integrands often arise 
in technical problems. In this paper, we assume that the integrand is 
an analytic function of the variable of integration u (in a suitable 
region) and tends to oscillate at a regular rate as«-> °° . The main 
concern here is how to deal with cases in which the rate of convergence 
is slow. 

A number of ways have been proposed to handle the slow con- 
vergence. One is to use integration formulas of Filon's type. In particu- 
lar, E. O. Tuck 2 has given a formula of this type suited to an infinite 
range of integration. Another method is to evaluate the contributions 
of positive and negative loops of the integrand and apply Euler's 
summation formula for slowly converging alternating series. Still 
another is to (i) write the integrand in terms of complex exponentials, 
(h) deform the path of integration so that it becomes a path (or 
portions of paths) of steepest descent in the complex M-plane, and then 
(Hi) integrate along straight-line segments that approximate the path. 
A closely related procedure is to tilt the path of integration, if possible, 
so that the integrand will decrease exponentially. 

155 



In this paper most of our attention will be devoted to still another 
method of improving the convergence, namely, that of subtracting 
from the integrand functions that behave like the leading terms in the 
asymptotic expansion of the integrand. This method is quite old. It is 
related to a method used by Kummer to deal with series, and its 
application to integrals is described by Krylov. 3 

Sections II through V are concerned with functions useful for 
subtraction when the limits of integration are and °o. Section VI 
deals with the evaluation of 

J(b,a) = e -6 / [exp (6w _1 sinw) — 1 — bu~ l sin w]e*° u dw, (1) 

J — ao 

which serves as a typical example of a slowly converging Fourier 
integral with limits ± «> . 

The integral (1) has been the subject of several papers because of its 
relation to the problem of calculating the power spectrum of a phase- 
modulated wave (see Prabhu and Rowe 4 and the references they 
give). In Section VI, results of trial calculations are presented that 
show that (1) can be efficiently evaluated by the trapezoidal rule after 
suitable functions have been subtracted from the integrand. 

Sections VII and VIII are concerned with the error introduced when 
a general Fourier integral with limits ± oo is evaluated by the trape- 
zoidal rule. 

II. PROPERTIES OF FUNCTIONS SUITABLE FOR SUBTRACTION- 
LIMITS AND oo 

Let the integral to be evaluated be 



1=1 f(u)du, 
Jo 



and let the leading term in the asymptotic expansion of f(u) for large 
u be u~" exp (iau), where a is real and v > 0. Quite often, in performing 
numerical integrations, this leading term is subtracted from the inte- 
grand only for values of u that exceed some large value U. We have a 
slightly different procedure in mind here. 

Here we want to use the trapezoidal rule methods described in Ref. 1. 
Since these methods require the integrand to be analytic, we assume 
that f(u) is analytic in a region containing the positive real axis. Our 
aim is to subtract an analytic function (which approaches u~* exp (iau) 
as u — * °o ) from f(u) over the entire interval (0, °° ). When this is done, 
the resulting integral has an integrand analytic over (0, » ) and can be 
evaluated by the trapezoidal rule methods described in Ref. 1 (see 
especially Section IX). 
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Functions suitable for subtraction when the limits are and oo 
should have the following properties : 

(i) As u — > oo , they should be asymptotic to u~" cos au or u~* sin au 
or, more generally, to 

w _ 'exp (iau), (2) 

where v > and a is real. 
(ii) They should be easy to compute and they should remain finite 

at u = 0. 
(tit) Their integrals between the limits (0, oo) should be easy to 

compute. 

In some cases, moderately simple functions can be used for sub- 
traction. If the leading term in the asymptotic expansion of the 
integrand is (sin au)/u, then we can subtract and add the integral 



/. 



u~ l sin au au = , ' . n (3) 

— ir/2, a < 0. 



Similarly, for (cos aw)/u we can use 

/•OO 

/ w _1 [cos au — exp ( — au)~\du = 0, (4) 

and, for (cos au)/u 2 , 

(cos au)du/(l + w 2 ) = £tt exp (— |o|). (5) 



/; 



The last integral (5) has the disadvantage that its integrand has the 
asymptotic expansion 

cos au cos au , ,„. 

-* ^-+-"' (6) 

and if we desire to subtract terms of both (u~ 2 ) and (u~*) from the 
original integrand, the term — (cos au)/u* in (6) has to be taken into 
account. In the applications we have in mind, the exp (—au) in (4) 
decreases so rapidly that it does not have to be taken into account 
in the way that — (cos aw)/u* does. 

III. GENERAL FUNCTIONS FOR SUBTRACTION— LIMITS AND co 

Here we present a set of functions that meets the "ease of compu- 
tation" requirements (ii) and (in) reasonably well. It is quite likely 
that there are better sets, but this one is the best I have been able to 
find. There are two formulas according to whether v is an integer or not. 
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The functions are the integrands in 

r ir— («— - e-^Ydu = r(1 , ~ 6) £ (-)-*( ? ) 

X[(n-fc)a + W+*" 1 , (7) 

X[(n - fc)a + *£Mfc[(n - k)a + fc/3], (8) 

where n is a positive integer and n, a, /3, and e are such that the integrals 
converge. In (7), e is not an integer, and (e) = 1, («)■ = e(e + 1) • • • 
(e + n — 1). In using (7) and (8), we set a = — ia/n, fi = \a\/n. 
In (7), n and e are chosen so that e < — 1 and n + e is equal to the 
exponent v in (2). The choice e < — 1 is made to make the integrand 
and its derivative in (7) approach as u — > 0. This tends to reduce the 
error in the numerical integration. Since o can be either positive or 
negative in a = — ia/n and /3 = | a \ /n, we take 

-tt/2 ^ arg [(n - k)a + *j8] ^ ir/2 

when calculating the logarithm and the (n + e — l)th power of 
l(n - k)a + kpj 

Equation (8) is obtained by letting e — > in (7), and (7) is based 
upon (page 243 of Ref. 5) 



•»—«-/;£ 






(9) 



where, in contrast to (7), e is restricted to < e < 1 to secure con- 
vergence. To obtain (7), expand [exp (—au) — exp (— /3w)] n by the 
binomial theorem. Then to each term exp []— (n — A;)aM — k(iu], add 
and subtract the series 

- E (-»)'[(» - *)« + W/li. 

Equation (7) then follows from (9) with t = [(n — k)a + A;/3]m and 
the relation 

L o (-)* ( I ) [(« - *)« + kfij - 0, (10) 

where Z = 0, 1, • • • , n — 1. Equation (10) expresses the fact that u n is the 
lowest power of u in the power series for [exp (—au) — exp (— /3w)] n . 
After (7) has been established for < e < 1, it can be extended by 
analytic continuation to the values of e needed in the application of (7). 
Other choices of a and /3 besides —ia/n and \a\/n can be made in 
(7) and (8). In some cases, it may be better to take a = — ia/n and 
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= a -\- b, where 6 is a positive constant at our disposal. For this 
choice of a and /3, [exp (—au) — exp (— /3w)] n becomes exp (iau) 
X [1 — exp( — 6w)] n and (n — k)a + fc/3 becomes kb — ia. By using 
these values in (7), equating imaginary parts, and letting e — > 1, we 
can obtain 

/ w -n-1 (l — e~ bu ) n sin au du 
Jo 

= Im -. t (-)""* ( ? ) ^b + to)*Z»(A;b + ia). 
n\k=o \k/ 

Here the integrand tends to (sin au)/u n+i as u — »*>. 

IV. EXAMPLE ILLUSTRATING APPLICATION OF EQ. (8) 

Consider the integral 

/ = r (1 + u 2 )-*e<"dt«, (11) 

which is equal to K (l) + tV[/ (l) - L (l)]/2, where K and / are 
Bessel functions and L is a Struve function (page 498 of Ref. 6). As 
u — > oo , the integrand approaches 

trV" - iw- 3 e ,u + • • • . 

Therefore, if we desire convergence like that of w _5 e iu , we apply (8) 
twice. First with n = 1, a = — i, = 1, and then with n = 3, 
a = -i/3,0 = *: 

,oo r / p iu _ p-u \ 1 / piult _ g-u/3 \ 3"] 

-^U?o(') ( - )3 " fcC(3 -" )( - i/3)+fc/3]2 

Xln[(3-ft)(-t/3)+*/3]. (12) 

The second line reduces to ztt/2 [which, incidentally, gives the com- 
bination of (3) and (4)]. Although the last sum appears complicated, 
it can be readily evaluated by a digital computer when the program 
is written using complex variables. 

The integrand in (12) now decreases as w _5 e ,u when u is large, and 
the integral can be evaluated by any suitable method of numerical 
integration. Some computations were made using the trapezoidal 
rule as described in Section IX of Ref. 1. In this method, the variable of 
integration is changed from u to x, where u = aln[l + exp (x/a)], 
and the integral in x is evaluated by the trapezoidal rule. A computa- 
tion using 52 points with h = Ax = 4.0 and a = 6 gave / = 0.42102 • • • 
+ iO. 87308- • •, which agrees with tabulated values. 
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This example brings to mind the question, when does one cease 
subtracting terms from the integrand and start integrating? That is, 
what is the trade-off between computing trapezoidal sums and those 
sums in (7) and (8)? There seems to be no well-defined answer, but 
experience indicates that there is usually no need to make the integrand 
decrease faster than the inverse fifth or sixth power of the variable of 
integration. 

Incidentally, the integral for / is quite readily evaluated by tilting 
the path of integration. Arbitrarily choosing a 45-degree tilt and setting 
u = (1 + i)t carries (11) into 

I = r eri{t)dt, 

f(t) = (l-M)(l-r-2^)-»e". (13) 

The contribution of the 45-degree arc at | u \ = » is zero by Jordan's 
lemma (page 115 of Ref. 5). The integral is now in the form discussed in 
Section VI of Ref. 1, and can be efficiently evaluated by the trapezoidal 
rule after the change of variable t = e v , v = x — e~ x . Computations 
made with this example suggest that when "tilting" can be used, it is 
preferable to "subtraction." However, tilting cannot always be used, 
an example being the integral J (6, a) defined by (1). 

V. EXAMPLE ILLUSTRATING THE USE OF EQ. (7) 

Consider the integral 

I = r (1 + u 2 )- 8/ V u dw, (14) 

the real part of which is related to the Bessel function K\iz{\). When 
u — > oo , the integrand becomes 

w - 6 /3 e .u _j_ 0(u~ nlz ). 

Therefore, to secure convergence of 0(u~ nlz ), we apply (7) once with 
n + e = 5/3. Since n must be a positive integer, possible choices of e 
are 2/3, -1/3, -4/3, •••. We also want e < - 1. We arbitrarily 
choose e = — 4/3 (and therefore n = 3) to make the integrand of (7) 
tend to as u ilz when u — > 0. Because n = 3 we take a = — i/ft and 
= 1/3. Then 



I = r du\_{\ + u 2 )- 6/6 e*' u - u- 5l3 (e iu ' a - e-»' 8 ) 8 ] 

M-4/3)(-l/3)(2/3)*V ; \k) 



X [(3-A0(-i73) + (fc/3)] 2 ' 3 . (15) 
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This integral can be evaluated by the trapezoidal rule in much the 
same way as was the integral in (12). 

VI. NUMERICAL EVALUATION OF THE INTEGRAL J(b,a) 

To simplify the discussion of Jib, a) defined by the integral (1), we 
assume that b is positive. The parameter a is real and J (b, a) is an 
even function of a. 

The integrand in J [b, a) can be expanded as exp (iau) times 

£ (6w -1 sin u) n /n\, 

and the term for n = 2 shows that the integral converges slowly. 
Inspection shows that the rate of convergence can be increased by 
writing J(b, a) as an integral from to a> and subtracting terms of 
the form vr n exp (icu) as described in Sections III and IV. However, 
it is more convenient to subtract terms of the form (w -1 sin u) " and 
use the known values of 

B„ = f {bu~ l sin u) n cos au du/n !, (16) 

where n ^ 1 is an integer and #„ is an even function of a. 

When n = 1, Bi is equal to br for \a\ < 1, to bir/2 for \a\ = 1, 
and to for \a\ > 1. When n ^ 2, B n is a continuous function of a 
which is when \a\ ^ n. When n ^ 2 and \a\ ^ n, 

'.-SS£ <->•(:)«■*- ^ (17) 

where M is the integer part of (n ± a)/2. The choice of the sign in =b 
is arbitrary, but it must be the same in M and in (17). The minus sign 
is the more convenient choice when a > 0. 

When we desire convergence at the rate of, say, 1/w 6 , we subtract 
and add e-"(B 2 + B 3 + B A + B & ) : 

J(b, a) = e~ b I exp (6w _1 sin u) — 1 

- E {bvr x sin u) n /rc!l e iau dw + e~ b (B 2 + £ 3 + B A + £ 5 ). (18) 

71=1 J 

Incidentally, J (6, a) is equal to e~ b times the infinite sum B 2 
+ £3 + • • • , a result useful for computation when b is small. 

The integrand in (18) is an analytic function of u, and the results 
of Ref. 1 suggest direct numerical integration by the trapezoidal rule. 
Several trial computations were made to obtain an idea of the perform- 
ance of this method of integration. In the actual work, the integral 
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Table I — Values of J(b,a) obtained by the trapezoidal rule 



b 


a 


h 


N 


Umax 


J(b, o) 


1 


1 


0.7 


12 


7.7 


0.413 5433 


1 


4 


0.4 


20 


7.6 


0.000 0428 


4 


1 


0.5 


36 


17.5 


1.341 1671 


4 


4 


0.4 


45 


17.6 


0.011 6253 


16 


1 


0.3 


39 


11.4 


0.997 3179 


16 


10 


0.225 


52 


11.5 


0.000 2046 


32 


1 


0.25 


19 


4.5 


0.736 6452 


32 


10 


0.175 


26 


4.4 


0.007 6251 



was rewritten as an integral from to °o , and exp (iau) was replaced 
by cos au. Typical results are shown in Table I. Here N is the number 
of terms in the trapezoidal sum, and h is the spacing between successive 
values of u. The trapezoidal sum was truncated at u = w mB x 
= (N — \)h. I believe the values shown for J(b, a) are accurate to 
the number of places shown. A short table of J(b, a) for a ^ 1.25 is 
given in Ref. 7, where J(b, a)e b is referred to as "Lewin's integral." 
When b is large, the work of Prabhu and Rowe 4 gives tight bounds on 
J (b, a) for all values of a. 

The truncation value w ma x was chosen to be the value of u beyond 
which the absolute value of the integrand remains less than 10 -8 . 
This makes the truncation error less than 10 -8 w max /5. It can be shown 
that w max has a maximum value of about 18 near 6 = 6. 

The error E in the computed value of J (b, a) resulting from the 
finite size of the spacing h, i.e., the "trapezoidal error," is of the order 
of J(b, 2irh~ 1 — a) when h is small. This follows from Section VII 
below. The trapezoidal errors observed in the construction of Table I 
agree well with the values of J(b, 2irhr l — a) estimated from bounds 
for J (b, a) given by Prabhu and Rowe. 4 E decreases rapidly as h 
decreases through the values shown in Table I. 

VII. TRAPEZOIDAL ERROR FOR FOURIER INTEGRALS 

Let the Fourier integral and the associated trapezoidal error be, 
respectively, 

1(a) = r Fivje^du, 

J —CO 



E = h £ F(nh) exp (ianh) -1(a). 
As in Ref. 1, Poisson's summation formula gives 

E = ( L + Z ) I(2Tkh~ 1 + a). 



(19) 
(20) 



(21) 
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When F(u) is an even function of u, I ( — a) is equal to 1(a), and 
(21) gives 

E = E UVirkh- 1 - a) + IVrkhr 1 + a)]. (22) 

It is interesting to note that, if 1(a) is known to be real and non- 
negative for all real values of a, then E is nonnegative and the values 
of 1(a) obtained by evaluating (19) by the trapezoidal rule are always 
greater than or equal to the true value of 1(a), truncation errors aside. 
This fact is of use in dealing with power spectra and probability 
densities. 

Now let F(u) be analytic in a strip in the w-plane containing the 
real w-axis. When F(u) is even and analytic, usually only the k = 1 
term in (22) is important if h is small. In addition, if a > 0, E and 
I(2irh~ l — a) are usually of the same order of magnitude. 

VIII. DISPLACING THE PATH OF INTEGRATION OF A FOURIER INTEGRAL 

In some cases, the numerical evaluation of the Fourier integral 1(a) 
denned by (19) can be facilitated by shifting the path of integration 
from the real w-axis to a parallel path that passes through u = ic, 
where c is some suitably chosen real constant. Changing the variable 
of integration in the integral (19) for 1(a) from u to x, where u = x + ic, 
and applying the trapezoidal rule to the integral in x gives the trape- 
zoidal error 

E' = h E F(nh + ic)e ianh - ac -1(a), 

which, by Poisson's summation formula, goes into 

E' = ( E + E ) ICZ-xkh-' + a)e i ^i h . (23) 

\A: = -« k=\/ 

When 1(a) is known to be real and nonnegative for all real values 
of a, then E' is nonnegative just as E is in (21), and the trapezoidal 
values of 1(a) obtained from the displaced path are greater than or 
possibly equal to the true value of 1(a), just as before. 

The integral in the expression (18) for J(b, a) is of the Fourier type 
(19), and some calculations were made to test the benefit obtained by 
shifting the path of integration. Along the path u = x + ic, the real 
part of the integrand is an even function of x, and the integral in (18) 
can be rewritten as an integral from x = to co of twice the real 
part. For b = 32 and a = 10, the value of c was chosen to be 0.87. This 
corresponds to the saddle point u = i0.87 of exp (bur 1 sin u + iau). 
A trapezoidal rule evaluation with h = 0.3 along the line u = x + t'0.87 
gave J(b, a) to the accuracy shown in Table I. The number of terms 
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required was N = 16 compared with iV = 26 when the integration 
was taken along the real w-axis (with h = 0.175). 
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